Competition between Normal Superfluidity and Larkin-Ovchinnikov Phases of 
Polarized Fermi Gases in Elongated Traps 
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By applying the recently proposed antisymmetric superfluid local density approximation (ASLDA) 
to strongly interacting polarized atomic gases at unitarity in very elongated traps, we find families of 
Larkin-Ovchinnikov (LO) type of solutions with prominent transversal oscillation of pairing poten- 
tial. These LO states coexist with a superfluid state having a smooth pairing potential. We suggest 
that the LO phase could be accessible experimentally by increasing adiabatically the trap aspect 
ratio. We show that the local asymmetry effects contained in ASLDA do not support a deformed 
superfluid core predicted by previous Bogoliubov-de Gennes treatments. 
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The s-wave pairing in two-component fermion systems 
is the most common mechanism inducing superconduc- 
tivity or superfluidity in the realm of condensed matter, 
nuclear physics, and cold atom gases [1|. The latter sys- 
tem, that was experimentally realized in recent years |2- 
y , has the great advantage that the pairing strength can 
be manipulated using Feshbach resonances to access the 
whole range of the BCS to BEC crossover. The competi- 
tion between s-wave pairing strength, trap deformation, 
population imbalance, and temperature could give rise 
to different superconducting phases and regions of phase 
separation. Recent experiments with two-component po- 
larized ultracold atomic gases in elongated traps open 
the possibility of the practical realization of long-sought 
exotic superfluid phases. Among them, the Fulde-Ferrer 
(FF) [7| state representing a condensate of Cooper pairs 
with a finite center-of-mass momentum and the Larkin- 
Ovchinnikov (LO) ^8] state with a spatially oscillating 
pairing potential have attracted much interest [9|, |lO| . 

The Rice experiment [5| on a polarized unitary gas 
has been at the center of an intense debate. They ob- 
served a spherical distortion of the superfluid core shape 
with respect to the elongated trap with an aspect ratio 
77^^50. Such a shape mismatch was not observed in the 
MIT experiment performed with a much less elongated 
trap (77~5) [J]. Recent experiments performed in traps of 
?7~22 [6| seem to confirm the MIT results. These seem- 
ingly conflicting experimental results have motivated sig- 



nificant theoretical efforts |llMl4| . 



Most of the theoretical approaches to trapped cold 
atoms have been based on the Local Density Approx- 
imation (LDA) or in the Bogoliubov-de Gennes mean- 
field approximation (BdG) with a contact interaction. 
These approximations, which are appropriate in the 
weak-coupling BCS region, are not expected to perform 
well in the strongly interacting region due to in-medium 



screening effects. However, most of the experiments are 
accomplished for strongly coupled dilute fermionic super- 
fluids in the unitary limit in which the s-wave scattering 
length approaches infinity {\as\ -^ oo). In this limit, the 
system exhibits universal behavior governed by the den- 
sities, making it ideally suited for a density functional 
theory (DFT) description [l5l-[l7|. DFT offers a way to 
incorporate experimental information, and the results of 
large-scale calculations like Monte Carlo [18|, Il9| within 
an accurate energy density able to describe strongly in- 
teracting systems due to the incorporation of many-body 
effects. 



Recently, SLDA [15|, |16|, a superfluid extension of 



DFT, was applied to symmetric two-component systems. 
Subsequently, it was generalized to polarized fermionic 
systems (ASLDA) [20| , finding strong evidence for the ex- 
istence of a stable LO phase in the thermodynamic limit. 
The LO phase was predicted to be stable in large regions 
of the parameter space of a polarized one-dimensional 
(ID) system '2l|, and very recently experimental evi- 
dence for a LO phase has been obtained in a spin mix- 
ture of ultracold ^Li atoms trapped in an array of ID 
tubes '221. Extremely elongated traps offer, therefore, 
a unique opportunity to explore the transition from 3D 
to ID polarized systems, and the competition between 
normal and exotic superfluid phases as a function of the 
aspect ratio. 

In this work we study cold polarized atoms in elon- 
gated traps at unitarity. The self-consistent BdG equa- 
tions of superfluid DFT are solved using the coordinate- 
space axial solver HFB-AX ^3, 24J, based on B-splines, 
which has been demonstrated to provide precise results 
at large deformations. In a previous work [23'| we showed 
the appearance of phase separation in deformed SLDA. 
Here, we compare SLDA results with those of the ASLDA 
formalism that has been tailored to describe imbalanced 



systems. Our focus is on phase separation effects, defor- 
mation of superfluid cores, and the appearance of super- 
fluid oscillating phases for different trap elongations. 

The grand-canonical energy density functional of 
ASLDA can be written as llBl 
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where the densities of spin- up p^(r) and spin-down pi{r) 
atoms, the pairing densities ^(r), and pairing gaps A(r) 
are: 
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K{r) =^f^u,{r)v*{r), A(r) =-ge//(r)K(r), 



(2) 



andp = pt + Pi- InEqs. (HE]), /« = 1/(1 + exp(^,/fcT)), 
Ui{r) and w,;(r) are eigenvectors of the BdG Hamilto- 
nian, and Ei is the corresponding eigenvalue. The BdG 
equations are solved self-consistently with the chemi- 
cal potentials Xa {<7—'[,i) being determined from the 
particle-number constraints N^y—TT^pcr) [12|, |25|, where 
iV-|- {N_i=N — iV-|-) is the number of spin- up (spin-down) 
fermions. The local polarization is given by x{r) = 
Pl{r)/p^{r) ^ 1 while the total polarization of the sys- 
tem is P={N^ — Ni)/N. The quantities aa-{x) and f3{x) 
that parametrize the local effective masses and normal 
interaction, respectively, were fitted to Monte Carlo data 
and experiment [20|. 

The minimization of the energy density ([1} leads to a 
BdG set of equations: 



h-f{r) — X^ 
A*(r) 



A(r) 
Mir) + Xi 



Ui{r) 
Vi{r) 



^E,, 



Ut{r) 
Vi{r) 



(3) 



where the single-particle Hamiltonian h„ = ~'Tm'^ ' 
(Vofo-) + Ua + Vext depends on the kinetic term, Hartree 
potential C/o-, and external harmonic trap potential Kxt- 
The ASLDA equations are solved with a finite energy 
cutoff Ec {\Ei\ ^ Ec) and a local regularized pairing 
interaction geffir) = 7[p^^'^(j") + A{r)j]~^ [16], where 
7 is the original pairing interaction, and A(r) is the 
average regulator for the spin-up and spin-down com- 
ponents. The SLDA formalism can be obtained from 
ASLDA assuming x{r) = 1, resulting in identical ef- 
fective masses and Hartree potentials for spin-up and 
spin-down species. The potential of the confining trap 
is K!xt(r) = ^{r"^ + z^ /rf), where 77 is the aspect ratio 
defining the trap elongation. The ASLDA equations are 
solved with hfb-AX in a discretized rectangular axial box 
(r, z). We work in trap units for which h ^ m = uj ^ 1. 
The energy cutoff was assumed to be Ec > SA-j- (in prac- 
tice, it was 35 in SLDA and 25 in ASLDA). 

We first study the occurrence of multiple superfluid 
solutions predicted in a BdG treatment of large-scale 
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FIG. 1: (Color online) SLDA results for a system of 1000 
atoms in a trap with aspect ratio r;=20 and polarization of 
P=0.37. (A) Energies of various self-consistent solutions of 
SLDA obtained by starting from oscillating initial conditions 
with different values of q (solid line). For each solution, the 
amplitude of the first oscillation in the pairing potential is 
shown (dashed line). (B) Pairing potentials for the three so- 
lutions with g=0 (I), g=1.05 (II), and q=2.2 (III). (C) Den- 
sity profiles p-f and pj, for the cases (I-III) along z (longitudi- 
nal direction; r = 0, solid line) and r\r (transverse direction; 
2 = 0, dotted line). Large deviations between longitudinal 
and transverse profiles are indicative of deformation eifects. 
If solid and dashed profiles are close, density distributions 
follow the geometry of the trap. See text for more details. 



trapped systems 11|. Similar to Ref. |ll| . we param- 



eterize an initial guess for the pairing potential along the 
z-axis as sin[g(z — Zc)]e~*-^~^'=-'/^, where Zc is the criti- 
cal distance at which oscillations develop. Figure [1] dis- 
plays SLDA results for Af=1000 atoms with polarization 
P=0.37 in a trap with 77=20. (It was advantageous to use 
SLDA in this systematic study of various self-consistent 
solutions as this approach is less computationally inten- 
sive than ASLDA.) We performed systematic calculations 
with different values of the wave number q flnding a mul- 
titude of self-consistent SLDA solutions. The panel A of 
Fig. [1] shows the energies of SLDA states as a function 
of q together with the amplitude of the flrst oscillation 
of the pairing potential. The ground-state (g.s.) solution 
has (7=0. In the energy curve, one can see two plateaus, 
the flrst one with 0.4<g<l.l containing nearly degener- 
ate solutions having large oscillation amplitudes and the 
second, with q > 1.5, containing nearly-degenerate solu- 
tions exhibiting low-amplitude oscillations. From them 
we select three typical solutions: (I) g.s. with 9=0, (II) 
excited state with g=1.05, and (III) excited state with 
(7=2.2. Panel B shows the corresponding pairing poten- 
tials along the z-axis and panel C displays the density 
proflles in longitudinal and transverse directions for the 
three states. Solution (II), with remarkable oscillations 



in the pairing potential and density profiles, is a LO state 
(see discussion below) . The higher-energy state (III) does 
not show appreciable oscillations in the density profiles. 
Unlike the solution I, solution II and particularly solu- 
tion III show the presence of the deformation effect in 
the minority component, similar to the one observed in 
the Rice experiment. 

The presence of nearly degenerate solutions at very 
large elongations of the trap should not be surpris- 
ing. Indeed, the trap excitation along the 2:-axis 
is Ai;^=fiwz/2=l/(277). In the example of Fig. [H 
AEz=0.025, which is '-0.5% of the Fermi energy. This 
means that the trap provides a quasi-continuum back- 
ground of states with smoothly increasing energy and 
wave number. We note that the wave number of state 



(II) qjj w qlo, where qlo 



is the LO 



wave number. All states belonging to the plateau around 
solution II have very similar values of chemical potentials, 
hence practically identical values of qlo- This suggests 
that the LO can be represented by a superposition of 
the DFT solutions with q around qlo UM- Indeed, it 
is anticipated that the inclusion of quantum fluctuations 
beyond DFT could mix nearly degenerate solutions, thus 
lowering the energy of the LO state with respect to the 
ground state. 

In Fig. [2] we compare density profiles in the longitudi- 
nal and transverse directions in SLDA and ASLDA for a 
g.s. configuration of 2000 polarized atoms with P=0.37 
in a moderately elongated trap of r]=10. The insets de- 
pict the pairing potential in the longitudinal direction. 
At this elongation, some low-amplitude oscillations in 
the pairing potential are present, but they are absent in 
the density profiles. By comparing the longitudinal and 
transverse density profiles, one can see a clear sign of de- 
formation of the core in the SLDA solutions as predicted 
by the BdG calculations of Refs. J12l.ll4[. However, when 
the local polarization effects are considered in ASLDA, 
the superconducting core follows the shape of the trap. 
Moreover, deformation effects in SLDA are washed out 
at temperatures around 20% of the Fermi temperature, 
see Fig. H 

Starting from the g.s. ASLDA solution at ?7=10, we 
evolve it gradually by increasing the aspect ratio. In 
this way we can reach an excited state of the system 
with q>Q. The resulting two-dimensional pairing poten- 
tials for ry=20, 30, and 40 are shown in Fig. [3l Figure [3] 
shows that the oscillations in the pairing potential, al- 
ready present for 77=20, grow dramatically with the as- 
pect ratio. The fact that the oscillations have transversal 
character (are radially aligned), also seen in large scale 
BdG calculations [11[ , gives another support for identify- 
ing this solution as a LO state. Moreover, the strength- 
ening of the oscillations at higher aspect ratios is consis- 
tent with the expectation of a LO phase in the quasi-lD 
limit [10|, |2l|, [22] . The numerical procedure for the solu- 
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FIG. 2: (Color online) Density profiles for the g.s. config- 
uration in SLDA (top) and ASLDA (bottom) in the longi- 
tudinal (solid lines) and transverse (dashed lines) directions. 
The pairing potentials A and the SLDA density profiles at 
temperature T=0.2 are shown in the insets. 




tion of ASLDA may simulate an experimental process in 
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FIG. 3: (Color online) Pairing potential A{r,z) evolved from 
the ASLDA g.s. solution of Fig.[2]to aspect ratios 77=20 (top), 
30 (middle), and 40 (bottom). The radial alignment of the 
nodes, characteristic of the LO state [ll[, is clearly seen. The 
LO wave numbers are g_Lo = 1.32, 1.18, and 1.13 for ^7=20, 30, 
and 40, respectively. Note the different scales in z and r. 



which a polarized atomic gas is formed in a moderately 
elongated trap and then the aspect ratio is increased adi- 
abatically to high elongations. 

In order to gain a deeper insight into the ASLDA for- 
malism, we show in Fig. |4] density profiles, Hartree po- 
tentials Ua, and effective masses a~^ in the longitudinal 
direction for the two atomic species in a trap of aspect 
ratio ?7=30 and P=0.37. All these quantities exhibit pro- 
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FIG. 4: (Color online) Densities, Hartree potentials, and effec- 
tive masses for an excited ASLDA state of 2000 atoms with 
P=0.37 in a trap with ri=30. The solid and dashed-dotted 
lines indicate the spin-up and spin-down components, respec- 
tively. 



nounccd oscillations consistent with the oscillations in 
the pairing potential of Fig. [3l Unlike in the SLDA 
formalism, these oscillations are very different for the 
two species. The Hartree potential, usually neglected in 
BdG calculations, has a significant infiuence on the self- 
consistent ASLDA solution. The rather shallow Hartree 
potential of the spin-down component does not favor the 
deformed-core solution, explaining the different density 
shapes seen in SLDA and ASLDA in Fig. [2] In these 
calculations, we used a large box with Zi„ax=90. We 
checked, however, that our results for profiles inside the 
box, are numerically stable with respect to small vari- 
ations in Zinax- Another interesting feature is that the 
effective masses, fitted to Monte Carlo calculations, are 
different at the trap boundaries. A recent experiment [6[ 
estimated the effective mass of the Fermi polaron to be 
1.17(10), while it was fitted to be 1.04 in Ref. ^. Our 
calculations with the new effective mass do not change 
appreciably from the presented ASLDA results. 

In summary, we have performed precise superfluid- 
DFT studies of spin-imbalanced Fermi gases in elongated 
traps at unitarity. Calculations using initial conditions 
with different pairing oscillations yield families of coex- 
isting pairing phases: the LDA solution characterized by 
a smooth pairing potential and LO type of solutions char- 
acterized by transversal oscillation of pairing potential in 
the area of phase separation. In a theory beyond BdG, 
such as the multireference DFT, these LO states can be 
viewed as superpositions of DFT solutions with different 
q- values. The deformed-core solutions predicted in SLDA 
for g.s. configurations of polarized systems in moderately 
elongated traps are absent in ASLDA. Consequently, the 
presence (or absence) of such states in experiment might 
provide clues about the effective interaction in strongly 
interacting polarized atomic gases. Our results suggest 
that LO could be studied experimentally by elongating 
the traps adiabatically, starting from a moderately elon- 
gated equilibrium state. 

Discussions with A. Bulgac and M.M. Forbes are grate- 



[1] R. Casalbuoni and G. NarduUi, Rev. Mod. Phys. 76, 263 

(2004); S. Giorgini, L.P. Pitaevskii, and S. Stringari, Rev. 

Mod. Phys. 80, 1215 (2008); D.J. Dean and M. Hjorth- 

Jensen, Rev. Mod. Phys. 75, 607 (2003). 
[2] W. Ketterle and M.W. Zwierlein, Riv. Nuovo Cimento. 

31, 247 (2008). 
[3] G.B. Partridge et ai, Science 311, 503 (2006). 
[4] Y. Shin et ai, Phys. Rev. Lett. 97, 030401 (2006). 
[5] G.B. Partridge et ai, Phys. Rev. Lett. 97, 190407 (2006). 
[6] S. Nascimbene et ai, Phys. Rev. Lett. 103, 170402 

(2009). 
[7] P. Fulde and R. A. FerreU, Phys. Rev. 135, A550 (1964). 
[8] A.I. Larkin and Y.N. Ovchinnikov, Sov. Phys. JETP 20, 

762 (1965). 
[9] T. Mizushima, K. Machida and M. Ichioka, Phys. Rev. 

Lett. 94, 060404 (2005); J. Kinnunen, L.M. Jensen, and 

P. Torma, Phys. Rev. Lett. 96, 110403 (2006); X.J. Liu, 

H. Hu, and P.D. Drummond, Phys. Rev. A 76, 043605 

(2007); M.R. Bakhtiari, M.J. Leskinen and P. Torma, 

Phys. Rev. Lett. 101, 120404 (2008); J.M. Edge and 

N.R. Cooper, Phys. Rev. Lett. 103, 065301 (2009); R.A. 

Molina, J. Dukelsky, and P. Schmitteckert, Phys. Rev. 

Lett. 102, 168901 (2009). 
[10] A. Bulgac and Michael McNeil Forbes, Phys. Rev. 

Lett. 101, 215301 (2008)^ 

[11] L.O. Baksmaty et al. , larXiv: 1003.4488 1 (2010) . 

[12] R. Sensarma et aL. larXiv:0706.1741ll (20071. 

[13] M. Ku, J. Braun, and A. Schwenk, Phys. Rev. Lett. 102, 

255301 (2009); T.N. De Silva and E.J. Mueller, Phys. 

Rev. Lett. 97, 070402 (2006); M. Haque and H.T.C. 

Stoof, Phys. Rev. Lett. 98, 260406 (2007); M.M. Parish 

and D.A. Huse, Phys. Rev. A 80, 063605 (2009). 

[14] M. Tezuka, Y. Yanase, and M. Ueda, a rXiv:0811. 16501 

(2010). 
[15] A. Bulgac and Y. Yu, Phys. Rev. Lett. 91, 190404 (2003). 
[16] A. Bulgac, Phys. Rev. A 76, 040502(R) (2007). 
[17] T. Papenbrock, Phys. Rev. A. 72, 041603(R) (2005). 
[18] J. Carlson et ai, Phys. Rev. Lett. 91, 050401 (2003). 
[19] G.G. Batrouni et al, Phys. Rev. Lett. 100, 116405 

(2008). 

[20] A. Bulgac and M.M. Forbes. [axXiv:0808.1436] (2008). 
[21] A.E. Feiguin and F. Heidrich-Meisner, Phys. Rev. B 76, 

220508(R) (2007); P. Kakashvih and C.J. Bolech, Phys. 

Rev. A 79, 041603(R) (2009); M. Tezuka and M. Ueda, 

(arXiv: 1002. 1433 (2010)^ 

[22] Yean-an Liao, et «;.. larXiv:0912.0092l (2009). 

[23] J.G. Pei, W. Nazarewicz, and M. Stoitsov, Eur. Phys. J. 

A 42, 595 (2009). 
[24] J.G. Pei et al, Phys. Rev. C 78, 064306 (2008). 
[25] G. Bertsch et al., Phys. Rev. A 79, 043602 (2009). 



